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We study the motion of small inertial particles in stratified turbulence. We derive a simplified 
model, valid within the Boussinesq approximation, for the dynamics of small particles in presence of 
a mean linear density profile. By means of extensive direct numerical simulations, we investigate the 
statistical distribution of particles as a function of the two dimensionless parameters of the problem. 

We find that vertical confinement of particles is mainly ruled by the degree of stratification, with 
a weak dependency on the particle properties. Conversely, small scale fractal clustering, typical of 
inertial particles in turbulence, depends on the particle relaxation time and is almost independent 
on the flow stratification. The implications of our findings for the formation of thin phytoplankton 
layers are discussed. 


Particles of density different from the surrounding fluid 
do not follow the motion of fluid particles, and generate 
inhomogeneous distributions even in incompressible flows 
[l[. This phenomenon is crucial in a variety of instances, 
from cloud formation in the atmosphere, to the dynamics 
of plankton in the ocean and lakes, to industrial appli¬ 
cations e.g. in reactors Q. Inhomogeneous distribution 
in turbulent flows is also of interest from a theoretical 
point of view. In recent years, analytical, numerical and 
experimental studies led to significant advances in the 


understanding of this process [3l-l8| . In an incompressible 
turbulent flow, non-inertial, fluid particles follow the flow 
streamlines and remain by definition homogeneously dis¬ 
tributed. In contrast, inertial particles with density dif¬ 
ferent from the fluid, are known to accumulate in regions 
of high vorticity (light particles) or high strain (heavy 


or mgn vorticity ( 
particles) iii, 


as a consequence of the accelerations 
induced by the flow. Recent analytical and numerical 
works have shown that gravity interacts with turbulent 
accelerations to increase clustering of inertial particles. 
Moreover, turbulence can increase the settling velocity 
with respect to still fluid, by pushing particles in regions 
of downward flow Elm. In presence of density fluctu¬ 
ations, gravity also affects the flow itself as in the case 
of stratified turbulence which finds many applications in 
natural and technological flows BP. One important 
example is ocean dynamics which is strongly affected by 
the presence of the pycnocline resulting from tempera¬ 
ture and salinity variations BB. 


Remarkably, very little is know about the distribution 
of inertial particles in stratified turbulence, in spite of 
its relevance for oceanic and other applications. Recent 
works have studied the effect of stratification on the clus¬ 
tering of heavy [l^ and light B particles and the ef¬ 
fect of a vertical confinement in homogeneous turbulence 
(B . In this Letter we investigate, by means of direct 
numerical simulations, the distribution of small buoyant 
particles transported in a turbulent stratified flow. We 
study both the large scale vertical distribution and the 


small scale clustering of particles as a function of the 
relevant parameters. Unexpectedly, we find that small 
scale (fractal) clustering is determined by the particles 
relaxation time and is almost independent on the degree 
of stratification. Conversely, the vertical confinement of 
particles is mostly controlled by the degree of fluid strat¬ 
ification, and weakly depends on the particle properties. 
The dependence of the confinement on Froude number is 
interpreted within the framework of stratified turbulence 
phenomenology. 



FIG. 1. Isopycnal surfaces 6 = z corresponding to the particle 
density po for two runs at Fr = 0.2 (upper plot) and Fr = 0.3 
(lower plot). The points represent the (fractal) distribution 
of particles with r = 0.01 in the stationary state. 


We consider a cubic box of size L of fluid linearly 
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(and stably) stratified in the direction 2 : of gravity g = 
( 0 , 0 , —g) with a constant mean density gradient dpjdz = 
— 7 . Within the Boussinesq approximation, the motion 
for the incompressible velocity field u(x, t) is ruled by 

O -I 

^+u-Vu = -—Vp + i/V^u-TV^ez + f (1) 
at po 


— +u-V9 = u-z + kS/^6 ( 2 ) 

together with V • u = 0. The scalar field which 

in the above equations has the dimension of a length, 
represents the deviations of the local density from the 
linear vertical profile, p = pq—^{z—9). v is the kinematic 
viscosity, k the density diffusivity and N = ( 73 /po)^^^ is 
the Brunt-Vaisala frequency, f represents an external 
mechanical forcing needed to sustain turbulence. In the 
inviscid, unforced limit (z/ = k = f = 0) equations (frill) 
conserve the total energy, sum of kinetic and potential 
contributions, E = i(|up) + where (...) denotes 

the average over the domain. 

The velocity v of a small inertial particle transported 
by the flow u generated by ([T]) evolves according to [l^ 


dt 




+ (1 - ^)g 


( 3 ) 


where /3 = 3p/(p + 2pp) is the density ratio {pp is the 
density of the particle of radius a) and Tp = a?/(2>i'j3) 
is the viscous relaxation time. We consider the limit of 
small particles (i.e. with small Tp for which one can ap¬ 
proximate ^ ~ ^ [ 1 ^) of density pp = po so that we 
can rewrite © as 


v = u-rp(l-/3)—-brp(l-/3)g (4) 

Consistently with the Boussinesq approximation, in the 
limit of small Tp, we can neglect the second term in the 
r.h.s. of (m and, since (1 — P)g ~ —^N^{z — 9)z, we ob¬ 
tain a simplified expression for the velocity of the floater 
whose position x evolves according to 

dx 

_=v = u--2-6Iz 5 

dt T 

where r = 3/{2N‘^Tp) represents the relaxation time 
of the particle onto the isopycnal surface of density 
p = Pq. This surface, defined implicitly by the relation 
2 = 9{x, y, 2 ), will be denoted as h{x, y), keeping in mind 
that in general it can be multivalued. We remark that 
more general models for the motion of floaters in strat¬ 
ified turbulence are possible, at the price of increased 
complexity and number of parameters na. 

Although the fluid velocity u is incompressible, the ve¬ 
locity field transporting the floaters is not since V • v = 
— (1 — d9/dz)/T which is in general nonzero. Formally, 


N 

Z/, K 

kf 

£l 

V 

Tp 

Re 

Fr 

T 

128 

5 

X 10“’^ 

1.0 

0.195 

0.028 

0.16 

430 

0.2 - 1.0 

0.01 

-10.0 

256 

4 

X 10“’^ 

1.0 

0.195 

0.024 

0.14 

540 

0.2 - 1.0 

0.01 

-10.0 


TABLE I. Parameters of the simulations. N resolution, z/ and 
K kinematic viscosity and diffusivity, kf forcing wavenumber, 
El energy input rate, 77 = Kolmogorov scale, Tp = 

Kolmogorov timescale. Re = UrmsL/v and Fr = 
/N. Re is computed at Fr = 1 at which the flow is 
almost isotropic, as Urms varies with Fr. The forcing scale 
is defined as L — tx jkf . All the simulations are performed at 
Schmidt number Sc = v j k = 1. 


this expression represents the rate of contraction of the 
phase space (here the configuration space) under the dy¬ 
namics. When it is negative we expect that trajectories 
of floaters will collapse on a (dynamical) fractal attrac¬ 
tor in the phase space. An example of the attractor is 
displayed in Fig. □ which shows that the large scale con¬ 
finement in the vertical direction coexists with a small 
scale clustering with fractal distribution on the isopycnal 
surface. 


We have integrated the Boussinesq equations dBUD in 
a domain of size L = 2tx with periodic boundary condi¬ 
tions by means of a fully parallel pseudo-spectral code 
at resolution up to = 256. Turbulence is generated 
by a (5-correlated in time isotropic forcing f which is ac¬ 
tive on a spherical shell of wavenumber around kf = 1 
and which pumps energy at the fixed rate e. These 
parameters define, together with the Brunt-Vaisala fre¬ 
quency, the Froude number Fr = {e^/^ky^)/N which 
measures the (inverse) stratification. We remark that, 
in stratified turbulence, the mean kinetic energy dissi¬ 
pation at small scales Sn is typically smaller than the 
input £ and depends on Fr, because a fraction of 
the energy input is converted into potential energy dur- 
ing the turbulent cascade and dissipated by diffusivity 
[20| . Another relevant parameter in stratified turbulence 
is the buoyancy Reynolds number Rcb = ej defined 
in terms of the ratio of the buoyancy (Ozmidov) scale 
to the dissipative scale 

as Rcb = in analogy to the usual Reynolds 

number Re = ■ These three numbers are not 

independent since Rcb = Fr^Re 2l| and Rcb = 1 dis¬ 
criminates between stratified-viscous flow {Rcb < 1 ) and 
stratified turbulence [Rcb > 1) [2l[. Our simulations 
are within the turbulent regime, as Rcb is in the range 
17 < Reb < 540. Together with the (HJH]), we integrated 
the equation (O for the particle motion for a set of 10 
classes of particles characterized by different values of r in 
the range 0.01 < t < 10.0. In presenting the results, this 
time will be made dimensionless with the Kolmogorov 
time Tp = (zz/e)^/^ by introducing a ’’Stokes number” 
St = t/t^. Table [J reports the most important parame¬ 
ters of the simulations. 

















3 


In Figure [2] we show vertical sections (at y = 0) of 
the isopycnal surface h (obtained from the solution of 
z = 0(x, y, z)) together with positions of the particles on 
the same sections, for different values of the parameters 
Ft and r. It is evident that the isopycnal surface h is 
almost flat for strong stratification and it becomes more 
bent (and multivalued) as Fr increases. Figure [2] shows 
also the effect of the relaxation time r on the particles. 
When t/tjj < I the particles are practically attached to 
the isopycnal surface, while their positions depart from 
the surface h by increasing r. 



FIG. 2. (Color online) Vertical [x, z) sections of isopycnal 
surfaces h (black points) for Fr = 0.3 (npper row) and Fr — 
1.0 (lower row) together with the positions of particles (green 
points) with St = 0.6 (left column) and St — 60.0 (right 
column). 


The probability density functions (PDF) of the verti¬ 
cal displacement of particles with respect to their equilib¬ 
rium position 2 = 0 in the absence of turbulence is shown 
in Fig. [3] for different values of Re, Fr and r. We found 
that in the wide range of parameters investigated these 
distributions are close to Gaussian (with some possible 
deviations in the tails). 

The standard deviations cr^ of these vertical distribu¬ 
tion of particles for different values of Fr and St are 
shown in Fig. U) We obtain a linear scaling of (7^ on Fr, 
with a coefficient which shows a (weak) dependency on 
St. For stronger stratification, Fr < 0.3, the standard 
deviation is almost independent on St, a feature which 
can be understood by looking at the plots in Fig. [21 It 
is indeed evident that for St not too large, the effect of 
the relaxation term in ([S]) is to allow the particle to de¬ 
tach from the level z = 9 and to remain “suspended” 
for a time of order r before feeling the vertical velocity 
towards the isopycnal surface. Because in stratified tur¬ 
bulence vertical velocity and vertical gradient of vertical 
velocity are suppressed 21| , there is no mechanism which 
ejects particles far from the surface. Therefore, as shown 



z/a^ 



z/o^ 


FIG. 3. Probability density functions of vertical positions of 
tracers for Re = 430, Fr = 0.3 (upper plot) and Re = 540?, 
Fr = 1.0 (lower plot). In both cases the PDF for different 
classes of floaters, corresponding to St = 0.06, 0.6, 6.0 and 
60.0 are shown. 



Fr 


FIG. 4. Standard deviations of the PDFs of particle ver¬ 
tical position as a function of Fr for different values of St 
at Re = 430. The solid black line represents the standard 
deviation cjh of the isopycnal surface. 
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by Fig. [2 the vertical region visited by particles reflects 
the vertical extension of the isopycnal surface and, there¬ 
fore, is independent on r. By increasing Fr, the vertical 
components of the velocity and of the velocity gradient 
increase and this produces a displacement of particles 
from the isopycnal surface when the relaxation time r is 
sufficiently large. 

In the limit of small 5't, the standard deviations of 
particles Oz collapse on the standard deviation Oh of the 
isopycnal surface. The linear dependence on Fr shown in 
Fig. |4]can be understood within the framework of strat¬ 
ified turbulence, as a manifestation of the presence of 
the so-called vertical shear layers 22| and the associated 
vertical correlation scale of velocity Physically this 
scale represents the vertical displacement for converting 
injected kinetic energy into potential energy and can be 
estimated simply as ~ 17 /N {U is a typical large scale 
velocity) and therefore one obtains Ly oc Fr, i.e. linear 
scaling as shown in Fig. |4l 





FIG. 5. Probability density functions of the variable ( = z —9 
for Re = 540, Fr = 1.0 and St = 1.4 (red line with square), 
St = 7.0 (blue line with circle) and St = 70 (pink line with 
triangle). The black line is a Gaussian. 

The deviation of the particles form the isopycnal sur¬ 
face can be investigated by looking at the statistics of 
the variable C = z — 9. The evolution of this quan¬ 
tity is obtained by the time evolution of the field 9 
along the trajectory of a floater moving with the ve¬ 
locity V and which gives (neglecting the diffusive term) 
d{z — 9)/dt = —{z — 9){1 — d9ldz)lT. In the absence of 
fluctuations {9 = 0) this equation would simply repre¬ 
sent the linear relaxation of particles towards the isopy¬ 
cnal layer z = 0 [l^. This is not achieved since the 
term d9ldz is fluctuating without a definite sign. Fig¬ 
ure [5] shows the normalized PDF of the variable C for 
three different values of r at Fr = 1.0. It is evident 
that the statistics is neither Gaussian nor scale invari¬ 
ant and the PDF develops large tails for small relaxation 
times. These large fluctuations are due to the folding of 


the isopycnal surface. Particles located in the neighbor¬ 
hood of a fold in which the isosurface is almost vertical 
(see examples in Figure [2|) are not restored horizontally 
to the close-by branch they just left, but rather displaced 
vertically by buoyancy toward the nearest branch above 
or below them. This mechanism, which is enhanced at 
large Fr, produces a sudden increase in the distance be¬ 
tween the particles and the isopycnal surface, and causes 
the development of large tails in the PDF of (. 



FIG. 6. Correlation dimension D 2 versus relaxation time St 
for different values of stratification parameter Fr. 

As already discussed, floating particles moving accord¬ 
ing to ([ 5 ]) are transported by a compressible velocity field 
and are therefore expected to relax on a (dynamical) frac¬ 
tal subset of the physical space, as shown in the exam¬ 
ples of Fig. „ In order to characterize this subset, and 
its dependence on the parameters, we have measured the 
correlation dimension D 2 of particle distribution, defined 
as the scaling exponent of the probability of finding two 
particles at distance less than r: F’dxi — X 2 I < r) oc r^^ 
as r —> 0 [ 2 ^ . The maximum value D 2 = 3 denotes 
uniformly distributed particles, while D 2 < 3 indicates 
fractal patchiness with smaller D 2 corresponding to more 
clustered distributions and increased probability of find¬ 
ing pairs of particles at close separation. Figure [6] shows 
D 2 as a function of St for different values of Fr. In the 
limit St ^ 1, for which v —> u in ([5]), floaters move as 
fluid particles in an incompressible velocity and therefore 
remain uniformly distributed in the volume with D 2 = 3. 
The fractal dimension is found to be monotonic in St and 
attains a minimum value D 2 — 1 for the smallest relax¬ 
ation time. This value indicates distributions of particles 
on quasi-one-dimensional structures (as shown qualita¬ 
tively in Fig.p) which is almost independent on Fr. This 
is a remarkable result, as one could expect that for strong 
stratification and small r, in which particles lie on the 
isopycnal surface which is almost flat, the fractal dimen¬ 
sion would be close to D 2 = 2. The fact that we find 
D 2 — 1 indicates that the dynamics on the surface is 
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dissipative and therefore the horizontal velocity field is 
compressible. 

The weak dependence of the correlation dimension on 
the stratification extends also to larger values of St and 
we find that, in general, D 2 is virtually independent on 
Fr for Fr > 0.5. This is in contrast with the (large scale) 
vertical confinement shown in Fig. [5] which displays a 
strong (linear) dependence on Fr and a weak dependence 
on St. 

In conclusion, we have investigated the dynamics of 
neutrally buoyant particles in a turbulent stratified flow. 
We have shown that the extension cr^ of the vertical layer 
in which the particles are confined depends on the char¬ 
acteristics of the flow and moderately on the size of the 
particles. On the contrary, the small scale patchiness in¬ 
side this layer, expressed by the correlation dimension, 
strongly depends on the particle size and only weakly on 
the stratification of the flow. 

One of the most remarkable examples of confinement 
of particles in the ocean is the formation of the so-called 
thin phytoplankton layers (TPL): aggregations of phy¬ 
toplankton and zooplankton at high concentration with 
thickness from centimeters to few meters, extending up to 
several kilometers horizontally and with timescale from 
hours to days [ 2 ^. 

Among the different mechanisms proposed for the for¬ 
mation of TPLs, buoyancy force in stratified flow is par¬ 
ticularly relevant for non-swimming species and aggre¬ 
gates, such as diatoms and marine snow, which are of¬ 
ten observed to accumulate in correspondence of strong 
stratification [ 2 ^. In order to test the applicability 
of our model in the aquatic ecosystems we provide a 
simple example using the typical values observed for 
diatom-dominated marine snow 2^. Assuming an en¬ 
ergy dissipation rate e ^ I0“® m^s“^ (r^ cs 10 s) and a 
Brunt-Vaisala frequency N ^ O.I s“^, aggregates of size 
a ~ 0.5 cm have a relaxation time r cs 18 s which corre¬ 
sponds to St = 1.8, inside the range in which we observe 
clustering at small scales (Fig.[6l) for all the values of Fr. 
Smaller aggregates or single cells of size a = O.I cm cor¬ 
respond to St ~ 45 and, according to our results, would 
distribute almost homogeneously (with a fractal dimen¬ 
sion close to 3) within the thin layer. 

Beside the applications to thin layers, our results are 
of general interest as they describe quantitatively and 
for the first time how the combination of turbulence and 


stratification generates both large scale confinement and 
small scale fractal patchiness in a suspension of buoyant 
particles. 

We acknowledge the European COST Action MP1305 
“Flowing Matter”. 
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